R Notebook

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(ggplot2)
library(forcats)
library(here)

here() starts at /Users/avrilwang/Desktop/Project-Plasmodium

library(deSolve)
library(crone)
library(optimParallel)

Loading required package: parallel

library(doParallel)

Loading required package: foreach
Loading required package: iterators

library(doRNG)

Loading required package: rngtools

library(arrow)

Attaching package: ‘arrow’

The following object is masked from ‘package:utils’:

    timestamp

library(stringr)
library(parallel)
library(ggpubr)

Notebook for plotting all of the figures for PloS Biology manuscript submission

Guidelines: taken from https://journals.plos.org/plosbiology/s/figures#loc-figure-file-requirements 1. format: tiff 2. max file size: 10 MB 3. text size: Arial, Times, or Symbol font only in 8-12 point 2. figure size: Width: 789 – 2250 pixels (at 300 dpi). Height maximum: 2625 pixels (at 300 dpi).

overall manuscript setup

1. Fig 1: model schematic
2. Fig 2. best single and co-infection infection cue [monocue] ——-> supplmentary table of parameters and fitness
3. Fig 3. dual cue optimization [monocue] ——-> supplementary table of flexible optimization
4. Fig 4. heterocue competition [heterocue]
5. Fig 5. partitioning single cue success and consequences of bottom+top down regulation [monocue]
6. Fig 6. consequences of adopting best cue [monocue]

#———————————–# best single and co-infection cue #———————————–# Figure displaying the reaction norms of best single and co-infection.

#——- optimal cue reaction norm ———–# # read data

# single infection dynamics, reaction norms, and rugs
si_dyn.df <- read_parquet(here("data/si_dyn/si_dyn_30.parquet")) 
si_rn.df <- read_parquet(here("data/si_dyn/si_rn.parquet"))
si_rug.df <- read_parquet(here("data/si_dyn/si_rug.parquet"))

# co-infection dynamics, reaction norms, and rugs
ci_dyn.df <- read_parquet(here("data/ci_dyn/ci_dyn.parquet"))
ci_rn.df <- read_parquet(here("data/ci_dyn/ci_rn.parquet"))
ci_rug.df <- read_parquet(here("data/ci_dyn/ci_rug.parquet"))

process data for reaction norm

# import labelling scheme
ez_label <- read.csv(here("data/ez_label.csv"))

# get si_label with ci cue range
si_ci_rug.df <- ci_rug.df %>% 
  mutate(label_si = case_when(
    label %in% c("I", "I1+I2") ~ "I",
    label %in% c("I log","I1+I2 log") ~ "I log",
    label %in% c("Ig", "Ig1+Ig2") ~ "Ig",
    label %in% c("Ig log") ~ "Ig log",
    label %in% c("sum", "I+Ig") ~ "I+Ig",
    label %in% c("sum log", "I+Ig log") ~ "I+Ig log",
    label == "R" ~ "R",
    label == "R log" ~ "R log",
    label %in% c("G", "G1+G2") ~ "G",
    label == "G log" ~ "G log"
  )) 

# get limit for si_rug
si_rug_lim.df <- si_rug.df %>% 
  filter(time <= 20) %>%
  group_by(label)%>% 
  summarise(min = min(value, na.rm = T)*0.9,
         max = max(value, na.rm = T)*1.1) %>% 
  select(label_si = label, min_si = min, max_si = max)

# filter to restriction conversion rate reaction norm range to cue ranges that appear in rug
## change to Inf/-inf to NA. Note that I am first joining with si rug lim to check which limit is larger, We will go with the cue range that has the largest span
ci_rug_lim.df <- si_ci_rug.df %>% 
  group_by(label) %>% 
  mutate(min = min(value, na.rm = T)*0.9,
         max = max(value, na.rm = T)*1.1) %>% 
  distinct(label, .keep_all = T) %>% 
  select(label, label_si, min, max)

rug_lim.final <- ci_rug_lim.df %>% left_join(si_rug_lim.df, by = "label_si") %>% 
  mutate(final_min = min(min, min_si),
         final_max = max(max, max_si))

# get second rug_lim.final for single infection
rug_lim.final2 <- rug_lim.final %>% 
  group_by(label_si) %>% 
  mutate(final_min = min(final_min, na.rm = T),
         final_max = max(final_max, na.rm = T))

# filter ci_rn by limit
ci_rn.df2 <- ci_rn.df %>% 
  left_join(rug_lim.final, by  = "label") %>% 
  group_by(label) %>% 
  filter(cue_range <= final_max & cue_range >= final_min)

match single infection rn with coinfection

# get ci label to si rug and filter by limit
si_rn.df2 <- left_join(si_rn.df, select(ez_label, label_si, label_ci), by = c("label" = "label_si")) %>% 
  left_join(rug_lim.final, by  = c("label_ci" = "label")) %>% 
  group_by(label_ci) %>% 
  filter(cue_range <= final_max & cue_range >= final_min) %>% 
  select(cue_range, cr, label_ci, label_si)

# get ci label to si rug
si_rug.df2 <- select(si_rug.df, value, label_si = label) 

plot reaction norm

# join with ezlabel
ci_rn.df3 <- ci_rn.df2 %>% left_join(ez_label, by = c("label" = "label_ci"))
si_rn.df3 <- si_rn.df2 %>% left_join(ez_label, by = "label_ci")
ci_rug.df3 <- ci_rug.df %>% left_join(ez_label, by = c("label" = "label_ci"))
si_rug.df3 <- si_rug.df2 %>% left_join(ez_label, by = "label_si")

# redo order of cues
ci_rn.df3$label <- factor(ci_rn.df3$label, 
                              levels = c("I", "I log",
                                         "I1+I2", "I1+I2 log",
                                         "Ig", "Ig log",
                                         "I+Ig", "I+Ig log",
                                         "sum", "sum log",
                                         "G", "G log",
                                         "G1+G2", "Ig1+Ig2",
                                         "R", "R log")) 

si_rn.df3$label_ci <- factor(si_rn.df3$label_ci, 
                              levels = c("I", "I log",
                                         "I1+I2", "I1+I2 log",
                                         "Ig", "Ig log",
                                         "I+Ig", "I+Ig log",
                                         "sum", "sum log",
                                         "G", "G log",
                                         "G1+G2", "Ig1+Ig2",
                                         "R", "R log"))

ci_rug.df3$label <- factor(ci_rug.df3$label, 
                              levels = c("I", "I log",
                                         "I1+I2", "I1+I2 log",
                                         "Ig", "Ig log",
                                         "I+Ig", "I+Ig log",
                                         "sum", "sum log",
                                         "G", "G log",
                                         "G1+G2", "Ig1+Ig2",
                                         "R", "R log"))

si_rug.df3$label_ci <- factor(si_rug.df3$label_ci, 
                              levels = c("I", "I log",
                                         "I1+I2", "I1+I2 log",
                                         "Ig", "Ig log",
                                         "I+Ig", "I+Ig log",
                                         "sum", "sum log",
                                         "G", "G log",
                                         "G1+G2", "Ig1+Ig2",
                                         "R", "R log"))

# plot
opt_cue_pl.A <- ggplot() +
  geom_line(data = ci_rn.df3, aes(x = cue_range, y = cr, color = "Co-infection")) +
  geom_line(data = si_rn.df3, aes(x = cue_range, y = cr, color = "Single infection"), linetype = "dashed") +
  geom_rug(data = ci_rug.df3, aes(x = value), color = "#4575b4", sides = "t") +
  geom_rug(data = si_rug.df3, aes(x = value), color = "#fc8d59", sides = "b") +
  facet_wrap(~ez_label, scales = "free_x", ncol = 2) +
  theme_bw() +
  labs(y = "Conversion rate", x = "Cue range", color = "Infection") +
  scale_x_continuous(labels = function(x) format(x, scientific = T),
                     guide = guide_axis(check.overlap = TRUE)) + 
                    theme(axis.text.x = element_text(size = 7))  +
  scale_color_manual(values=c( "#4575b4", "#fc8d59"))

ggsave(here("figures/plos-bio/reaction_norm.png"), width = 7, height = 10)

#———- plot best fitness ————–# # import in data

# single infection dynamics
si_dyn.df <- read_parquet(here("data/si_dyn/si_dyn_30.parquet")) 

# co-infection dynamics
ci_dyn.df <- read_parquet(here("data/ci_dyn/ci_dyn.parquet"))

ez_label <- read.csv(here("data/ez_label.csv"))

process for final 20 days fitness

# get single infection maximum tau_cum for 20 days
si_fitness.df <- si_dyn.df %>% 
  filter(variable == "tau_cum" & time == 20)
# get single infection maximum tau_cum for 20 days
si_fitness.df <- si_dyn.df %>% 
  filter(variable == "tau_cum" & time == 20)
# get co-infection maximum tau_cum for 20 days
ci_fitness.df <- ci_dyn.df %>% 
  filter(variable == "tau_cum1" & time == 20)

Error in filter(., variable == "tau_cum1" & time == 20) : 
  object 'ci_dyn.df' not found

plot

opt_cue_pl.B <- ggplot() +
  geom_point(data = si_ci_fitness.df, aes(x = fitness_si, y = fitness_ci, color = ez_label_si, shape = ez_label_si), size = 2.5) +
  ggrepel::geom_label_repel(data = si_ci_fitness.df, aes(label = ez_label, x = fitness_si, y = fitness_ci),
                            fill = "white",xlim = c(-Inf, Inf), ylim = c(NA, NA)) +
  labs(x = "Optimal single-infection fitness", y = "Optimal Co-infection fitness (per strain)",
       color = "Single infection cue", shape = "Single infection cue") +
  scale_shape_manual(values = 15:24) +
  lims(x = c(8, 10.5)) +
  geom_vline(xintercept = 9.2, linetype = "dashed", alpha = 0.5) + # si good cue cutoff
  geom_hline(yintercept = 2.25, linetype = "dashed", alpha = 0.5) + # ci good cue cutoff
  theme_bw() +
  coord_cartesian(clip = "off") +
  theme(plot.margin=margin(r=0))

#————– timeseries conversion rate —————–# # process info for single infection

plot single infection conversion rate heatmap

# plot poor performing
si_cr.pl1 <- ggplot() +
  geom_tile(data = si_cr.poor, aes(x = time, y = forcats::fct_reorder(ez_label_si, fitness_si), fill = value)) +
  labs(x = "Time (days)", y = "Low performing\nsingle infection cues", fill = "Conversion rate") +
  scale_fill_viridis_c(limits = c(0, 1)) +
  xlim(1, 20) +
  theme_bw()

# plot high perfomring
si_cr.pl2 <- ggplot() +
  geom_tile(data = si_cr.high, aes(x = time, y = forcats::fct_reorder(ez_label_si, fitness_si), fill = value)) +
  labs(x = "", y = "High performing\nsingle infection cues", fill = "Conversion rate") +
  scale_fill_viridis_c(limits = c(0, 1)) +
  xlim(1, 20) +
  theme_bw()

save figure for poster

ggarrange(si_cr.pl2, si_cr.pl1, common.legend = T, ncol = 1, nrow = 2)

Warning: Removed 6000 rows containing missing values (geom_tile).
Warning: Removed 6000 rows containing missing values (geom_tile).
Warning: Removed 10000 rows containing missing values (geom_tile).

co-infection conversion rate heatmap

plot co-infeciton convesion rate heatmap

plot

# plot poor performing
ci_cr.pl1 <- ggplot() +
  geom_tile(data = ci_cr.poor, aes(x = time, y = forcats::fct_reorder(ez_label_si, fitness_ci), fill = value)) +
  labs(x = "Time (days)", y = "Low performing\nco-infection cues", fill = "Conversion rate") +
  scale_fill_viridis_c(limits = c(0, 1)) +
  xlim(1, 20) +
  theme_bw()

# plot high perfomring
ci_cr.pl2 <- ggplot() +
  geom_tile(data = ci_cr.high, aes(x = time, y = forcats::fct_reorder(ez_label_si, fitness_ci), fill = value)) +
  labs(x = "", y = "High performing\nco-infection cues", fill = "Conversion rate") +
  scale_fill_viridis_c(limits = c(0, 1)) +
  xlim(1, 20) +
  theme_bw()

#——— assemble final figure ————–#

tiff(here("figures/plos-bio/opt_cue.tiff"), units="px", width = 2250, height = 1500, res=300, compression = "lzw", type="quartz")
ggarrange(opt_cue_pl.A, opt_cue_pl.BC, ncol = 2, nrow = 1, labels = c("A", ""), widths = c(0.75, 1))

Warning: ggrepel: 7 unlabeled data points (too many overlaps). Consider increasing max.overlaps

dev.off()

null device 
          1 

#——————————————–# # dual cue optimization figure #——————————————–#

source(here("functions/chabaudi_si_clean_high.R"))
source(here("functions/chabaudi_si_clean.R"))
source(here("functions/par_to_hm_te.R"))

#———- plotting fitness of dual vs single cue opt ———# # import in previous data

# dual cue fitness
dual_fitness.df <- read.csv(here("data/dual_cue_opt4/dual_cue_fitness_20.csv"))
## make label and filter out very low fitness
dual_fitness.df <- dual_fitness.df %>% 
  mutate(temp_label = gsub("log", "log10", label),
         temp_label_b = gsub("log", "log10", label_b),
         label_final = paste0(temp_label, "+", temp_label_b)) %>% 
  filter(value > 2)

# get single cue fitness
si_dyn.df <- read_parquet(here("data/si_dyn/si_dyn_30.parquet")) 
si_fitness.df <- si_dyn.df %>% 
  filter(variable == "tau_cum" & time == 20)

# join si and dual cue
dual_si_fitness.df <- dual_fitness.df %>% 
  left_join(select(si_fitness.df, id, si_fitness = value), by = "id") %>% 
  left_join(select(si_fitness.df, id_b = id, si_fitness_b = value), by = "id_b") %>% 
  mutate(si_fitness_max = ifelse(si_fitness > si_fitness_b, si_fitness, si_fitness_b),
         dual_label = gsub("log", "log10", paste(label, "+", label_b)))

plot

#———– time series conversion rate ————-# # dynamics simulation of high parameter cues (these serve as reference points)

# best dual cue combo
dual.cr <- chabaudi_si_clean(
  parameters_cr = c(4.446192033,    10.97518275,    1.38762817, 23.3059254, -3.452052371,   -18.0070692,    39.66614226,    -3.545193141,   18.78350799),
  immunity = "tsukushi",
  parameters = parameters_tsukushi,
  time_range = seq(0, 20, by = 1e-3),
  cue_range =  seq(6, 7, by = 1/500),
  cue_range_b = seq(0, log10(6*(10^6)), by = (log10(6*(10^6)))/500),
  cue = "R",
  cue_b = "I",
  log_cue = "log10",
  log_cue_b = "log10",
  solver = "vode",
  dyn = T
)

# when time is used as a cue (high parameter)
time.cr <- chabaudi_si_clean_high(
  parameters_cr = c(9.154314,  -7.570829, -22.506638 ,  3.382405 ,-13.453519 ,-17.011485  , 3.678181, -12.851895 ,-26.115158),
  immunity = "tsukushi",
  parameters = parameters_tsukushi,
  time_range = seq(0, 20, by = 1e-3),
  cue_range =  seq(0, 20, by = 1e-3),
  cue = "t",
  solver = "vode",
  dyn = T)

# when asexual iRBC is used as a cue (high flexibility)
I_high.cr <- chabaudi_si_clean_high(
  parameters_cr = c(1.296675,  3.544034 , 4.907484,  2.174249, -3.238309 ,-5.181614 ,-1.645072 , 1.834302 , 1.581011),
  immunity = "tsukushi",
  parameters = parameters_tsukushi,
  time_range = seq(0, 20, by = 1e-3),
  cue_range =  seq(0, log10(6*(10^6)), by = (log10(6*(10^6)))/5000),
  cue = "I",
  log_cue = "log10",
  solver = "vode",
  dyn = T)

# when asexual iRBC is used as a cue (normal flexibility)
I.cr <- chabaudi_si_clean(
  parameters_cr = c(5.463558,   2.383948,   -17.757281, 4.571835),
  immunity = "tsukushi",
  parameters = parameters_tsukushi,
  time_range = seq(0, 20, by = 1e-3),
  cue_range =  seq(0, log10(6*(10^6)), by = (log10(6*(10^6)))/5000),
  cue = "I",
  log_cue = "log10",
  solver = "vode",
  dyn = T)

# process 
I_high.cr2 <- I_high.cr %>% filter(variable == "cr") %>% mutate(label_new = "I log10 flexible") %>% select(-variable)

I.cr2 <- I.cr %>% filter(variable == "cr") %>% mutate(label_new = "I log10") %>% select(-variable)

time_high.cr2 <- time.cr %>% filter(variable == "cr") %>% mutate(label_new = "Time flexible") %>% select(-variable)

dual.cr2 <- dual.cr %>% filter(variable == "cr") %>% mutate(label_new = "R log10 + I log10") %>% select(-variable)

# combine
dual_cr.df <- rbind(I_high.cr2, I.cr2, time_high.cr2, dual.cr2)

plot

ggplot() +
  geom_line(data = dual_cr.df, aes(color = label_new, x = time, y = value), size = 1) +
  geom_point(data = dual_cr.df %>% filter(time%%1 == 0), aes(color = label_new, x = time, y = value, shape = label_new), size = 2) +
  labs(x = "Time (days)", y = "Conversion rate", color = "Cue", shape = "Cue") +
  xlim(0, 20) +
  scale_color_manual(values = c("#4575b4", "#91bfdb","#fc8d59","#fdcb44")) +
  theme_bw() +
  theme(legend.position="bottom") +
  guides(color = guide_legend(nrow = 2, byrow = TRUE))

Warning: Removed 3 row(s) containing missing values (geom_path).
Warning: Removed 3 rows containing missing values (geom_point).

#———— reaction norm heatmap of R log10 + I log10 ————# # process data

# make heatmap df
R_I.hm <- par_to_hm_te(par = c(4.446192033, 10.97518275,    1.38762817, 23.3059254, -3.452052371,   -18.0070692,    39.66614226,    -3.545193141,   18.78350799),
             cue_range = seq(6, 7, length.out = 500),
             cue_range_b = seq(0,   6.77815125, length.out = 500))

Error in par_to_hm_te(par = c(4.446192033, 10.97518275, 1.38762817, 23.3059254,  : 
  could not find function "par_to_hm_te"

plot

ggplot() +
  geom_raster(data = R_I.hm, aes(x = cue_range_b, y = cue_range, fill = cr)) +
  scale_fill_viridis_c() +
  geom_path(data = R_I.dyn, aes(x = log_I, y = log_R), color = "white", arrow = arrow(angle = 30, length = unit(0.1, "inches"))) +
  geom_point(data = R_I.dyn %>% filter(row_number() %% 1000 == 1 & time <= 20), aes(x = log_I, y = log_R), color = "white") +
  xlim(0.99*min(hablar::s(R_I.dyn$log_I), na.rm = T), 1.01* max(hablar::s(R_I.dyn$log_I), na.rm = T)) +
  ylim(0.99*min(hablar::s(R_I.dyn$log_R), na.rm = T),1.01* max(hablar::s(R_I.dyn$log_R), na.rm = T)) +
  labs(y = "RBC log10", x = "Asexual iRBC log10", fill = "Conversion rate") +
  theme_dark()

Warning: Removed 219202 rows containing missing values (geom_raster).

save figure for poster

ggarrange(dual_cr.pl2, dual_rn.pl2, align = "h", widths = c(1.25, 1))

Warning: Removed 3 row(s) containing missing values (geom_path).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 219202 rows containing missing values (geom_raster).

#——– assemble final figure ————-#

# assemble panel A
dual.pl <- ggarrange(dual_si_fitness.pl, dual_pl.BC, ncol = 2, labels = c("A", ""), width = c(1, 0.75))

Warning in as_grob.default(plot) :
  Cannot convert object of class numeric into a grob.

#——————————————-# # heterocue adoption #——————————————-#

#——- static competition —————# # calculate fitness difference for 20 days

import and process data

min(static.df4$fitness_difference)

[1] -0.922165

plot

static.pl <- ggarrange(static.pl1, static.pl2, align = "hv", widths = c(1, 0.2))

Warning: Graphs cannot be vertically aligned unless the axis parameter is set. Placing graphs unaligned.

ggsave(here("figures/plos-bio/static_heatmap.tiff"), width = 10, height = 6)

#—— invasion matrix ——————–# # import in data (already 20 days )

invade.df <- read.csv(here("data/ci_invasion.csv"))

process data for invasion matrix

invade.mat4 <- rbind(invade.mat3, 
                        invade.mat3 %>% 
                          rename(V2_label = V1_label, V1_label = V2_label) %>% 
                          select(V1_label, V2_label, invasion)) %>%
  mutate(across(ends_with("label"),
                ~factor(.x, levels = levels))) %>%
  filter(as.integer(V1_label) < as.integer(V2_label)) %>%
  mutate(V2_label = forcats::fct_rev(V2_label)) 

Adding missing grouping variables: `id_alt`

plot invasion matrix

create summary bar chart

# create a stacked barchart for summary
## filter out na
invade.matalt <- invade.mat3 %>% na.exclude()

# get frquency from both sides. Note when grouping for V2, from the perspective of cue 2, scenarrio when strain 2 invade = strain 1 invade
invade.matalt1 <- invade.matalt %>% group_by(V1_label, invasion) %>% 
  summarize(frequency_1 = n())

invade.matalt2 <- invade.matalt %>%
  mutate(invasion_alt = case_when(
    invasion == "Only strain 1 invasion" ~ "Only strain 2 invasion",
    invasion == "Only strain 2 invasion" ~ "Only strain 1 invasion",
    invasion == "Mutual invasion" ~ "Mutual invasion",
    invasion == "Mutual non-invasion" ~ "Mutual non-invasion"
  )) %>% 
  group_by(V2_label, invasion_alt) %>% 
  summarize(frequency_2 = n())     

# full join and sum. has confirmed all of them add up to 14 
invade.matalt3 <- full_join(invade.matalt1, invade.matalt2, by = c("V1_label" = "V2_label", "invasion" = "invasion_alt"))

invade.matalt3[is.na(invade.matalt3)] <- 0
invade.matalt4 <- invade.matalt3 %>% 
  mutate(freq = frequency_1 + frequency_2) %>% 
  mutate(temp = case_when(
    invasion == "Only strain 1 invasion" ~ freq
  )) %>% 
  group_by(V1_label) %>% 
  mutate(invade_1_freq = max(temp, na.rm = T))
  
invasion.pl2 <- ggplot() +
  geom_bar(data = invade.matalt4, aes(x = freq, y = reorder(V1_label, invade_1_freq), fill = forcats::fct_rev(factor(invasion, levels = c("Only strain 1 invasion", "Only strain 2 invasion", "Mutual invasion", "Mutual non-invasion")))), stat = "identity") +
  labs(x = "Frequency", fill = "Invasion", y = "Strain 1 cue") +
  theme_bw()  +
  scale_fill_manual(values = c("Only strain 1 invasion" = "#4575b4",
                               "Only strain 2 invasion" = "#fc8d59",
                               "Mutual invasion" = "#fee090")) +
  theme(text = element_text(size = 15))

plot together

ggsave(here("figures/plos-bio/invasion_heatmap.tiff"), width = 10)

Saving 10 x 7 in image

#—— effect cue perception ——-#

static

logging

# get non-logged pairings
static_nolog <- static.df2 %>% 
  mutate(cue_1 = trimws(gsub("\\ .*", "", label_ci_1)),
         log_1 = case_when(
           str_detect(label_ci_1, "log") ~ "log",
           str_detect(label_ci_1, "log", negate = T) ~ "none")) %>% 
  filter(log_1 == "none")


static_log <- static.df2 %>% 
  mutate(cue_1 = trimws(gsub("\\ .*", "", label_ci_1)),
         log_1 = case_when(
           str_detect(label_ci_1, "log") ~ "log",
           str_detect(label_ci_1, "log", negate = T) ~ "none")) %>% 
  filter(log_1 == "log")

static_log.df <- left_join(
  select(static_nolog, cue_1, label_ci_2, log_1, None = fitness_diff),
  select(static_log, cue_1, label_ci_2, log_1, Log = fitness_diff),
  by = c("cue_1", "label_ci_2")) %>% 
  filter(!is.na(None) & !is.na(Log)) %>% 
  mutate(classification = ifelse(Log > None, "Logged better", "Not logged better"))

Error in `stop_subscript()`:
! Can't subset columns that don't exist.
x Column `fitness_diff` doesn't exist.
Backtrace:
  1. ... %>% ...
  6. dplyr:::select.data.frame(...)
  7. tidyselect::eval_select(expr(c(...)), .data)
  8. tidyselect:::eval_select_impl(...)
 16. tidyselect:::vars_select_eval(...)
     ...
 23. tidyselect:::chr_as_locations(x, vars)
 24. vctrs::vec_as_location(x, n = length(vars), names = vars)
 25. vctrs `<fn>`()
 26. vctrs:::stop_subscript_oob(...)
 27. vctrs:::stop_subscript(...)

combined

static_comb.df <- left_join(
  select(static_nocomb, cue_1, label_ci_2, log_1, Self = fitness_diff),
  select(static_comb, cue_1, label_ci_2, log_1, Total = fitness_diff),
  by = c("cue_1", "log_1", "label_ci_2")) %>% 
  filter(!is.na(Total) & !is.na(Self)) %>% 
  mutate(classification = ifelse(Total > Self, "Total better", "Self better"))

Error in `stop_subscript()`:
! Can't subset columns that don't exist.
x Column `fitness_diff` doesn't exist.
Backtrace:
  1. ... %>% ...
  6. dplyr:::select.data.frame(...)
  7. tidyselect::eval_select(expr(c(...)), .data)
  8. tidyselect:::eval_select_impl(...)
 16. tidyselect:::vars_select_eval(...)
     ...
 23. tidyselect:::chr_as_locations(x, vars)
 24. vctrs::vec_as_location(x, n = length(vars), names = vars)
 25. vctrs `<fn>`()
 26. vctrs:::stop_subscript_oob(...)
 27. vctrs:::stop_subscript(...)

plot

static_log.pl <- ggpaired(static_log.df, cond1 = "None", cond2 = "Log", line.color = "classification", alpha = 0.5) +
  labs(x = "Cue perception", y = "Fitness difference", color = "Effect") +
  scale_color_manual(values = c("Not logged better" = "#fc8d59", "Logged better" = "#4575b4"))

static_comb.pl <- ggpaired(static_comb.df, cond1 = "Total", cond2 = "Self", line.color = "classification", alpha = 0.5) +
  labs(x = "Cue perception", y = "Fitness difference", color = "Effect") +
  scale_color_manual(values = c("Total better" = "#fc8d59", "Self better" = "#4575b4"))

static_cue.pl <- ggarrange(static_log.pl, static_comb.pl, ncol = 2, align = "h")
ggsave(here("figures/plos-bio/static_cue_perception.tiff"), width = 7.5, height = 4)

invasion

proces data

# join invade df with label because I am lazy
invade.df2 <- invade.df %>% 
  left_join(select(ez_label, id_ci, label_ci_1 = label_ci), by = c("mut_id" = "id_ci"))

log

# get non-logged pairings
invade_nolog <- invade.df2 %>% 
  mutate(cue_1 = trimws(gsub("\\ .*", "", label_ci_1)),
         log_1 = case_when(
           str_detect(label_ci_1, "log") ~ "log",
           str_detect(label_ci_1, "log", negate = T) ~ "none")) %>% 
  filter(log_1 == "none")


invade_log <- invade.df2 %>% 
  mutate(cue_1 = trimws(gsub("\\ .*", "", label_ci_1)),
         log_1 = case_when(
           str_detect(label_ci_1, "log") ~ "log",
           str_detect(label_ci_1, "log", negate = T) ~ "none")) %>% 
  filter(log_1 == "log")

invade_log.df <- left_join(
  select(invade_nolog, cue_1, res_id, log_1, None = fitness),
  select(invade_log, cue_1, res_id, log_1, Log = fitness),
  by = c("cue_1", "res_id")) %>% 
  filter(!is.na(None) & !is.na(Log)) %>% 
  mutate(classification = ifelse(Log > None, "Logged better", "Not logged better"))

combined

invade_nocomb <- invade.df2 %>% 
  mutate(cue_1 = ifelse(
    str_detect(label_ci_1, "sum"), "I+Ig",
    trimws(gsub("+1.*|\\ log", "", label_ci_1))
  ),
  log_1 = case_when(
           str_detect(label_ci_1, "log") ~ "log",
           str_detect(label_ci_1, "log", negate = T) ~ "none"),
    comb_1 = case_when(
    str_detect(label_ci_1, "1+|sum") ~ "comb",
    str_detect(label_ci_1, "1+|sum", negate = T) ~ "none" 
  )) %>% 
  filter(comb_1 == "none")

invade_comb <- invade.df2 %>% 
  mutate(cue_1 = ifelse(
    str_detect(label_ci_1, "sum"), "I+Ig",
    trimws(gsub("+1.*|\\ log", "", label_ci_1))
  ),
  log_1 = case_when(
           str_detect(label_ci_1, "log") ~ "log",
           str_detect(label_ci_1, "log", negate = T) ~ "none"),
    comb_1 = case_when(
    str_detect(label_ci_1, "1+|sum") ~ "comb",
    str_detect(label_ci_1, "1+|sum", negate = T) ~ "none" 
  )) %>% 
  filter(comb_1 == "comb")
  
invade_comb.df <- left_join(
  select(invade_nocomb, cue_1, res_id, log_1, Self = fitness),
  select(invade_comb, cue_1, res_id, log_1, Total = fitness),
  by = c("cue_1", "log_1", "res_id")) %>% 
  filter(!is.na(Total) & !is.na(Self)) %>% 
  mutate(classification = ifelse(Total > Self, "Total better", "Self better"))

plot

invade_cue.pl <- ggarrange(invade_log.pl, invade_comb.pl, align = "h", ncol = 2)
ggsave(here("figures/plos-bio/invasion_cue_perception.tiff"))

Saving 7 x 7 in image

#——————————————-# # Partitioning best cue #——————————————-# #——- single infection ———–# # redo some optimization (lower fitness in no R than default)

source(here("functions/chabaudi_si_clean_R.R"))
source(here("functions/chabaudi_si_clean_N.R"))
# I none
cl <- makeCluster(detectCores()); setDefaultCluster(cl = cl)
I_no_R <- optimParallel(
    par = rep(0.5,4), # start at 0.5x4
    fn = chabaudi_si_clean_R, 
    control = list(trace = 6, fnscale = -1),
    immunity = "tsukushi",
    parameters = parameters_tsukushi,
    time_range = seq(0, 20, by = 1e-3),
    cue_range =  seq(0, 6*(10^6), by = (6*(10^6))/5000),
    cue = "I",
    log_cue = "none",
    solver = "vode")
stopCluster(cl)
# 0.144021 -43.1046 2030.27 -524.686 
# 8.69589

import and process data

# import in data
si_partition.ls <- list.files(path = here("data/partition/si/"), pattern = "*.csv", full.names = T)
si_partition.ls <- lapply(si_partition.ls, read.csv)

Warning in read.table(file = file, header = header, sep = sep, quote = quote,  :
  incomplete final line found by readTableHeader on '/Users/avrilwang/Desktop/Project-Plasmodium/data/partition/si//I_none_partition.csv'

si_partition.df <- do.call(rbind, si_partition.ls)

# combine with si fitness (default)
si_partition.df <- si_partition.df %>% left_join(select(si_fitness.df, id, fitness = value), by = "id")

Error in select(si_fitness.df, id, fitness = value) : 
  object 'si_fitness.df' not found

plot

# raw fitness
si_partition.pl1 <- ggplot() +
  geom_vpline(data = si_partition.df2, aes(y = fct_reorder(category, mean), x = mean, group = category, color = category), show.legend = F, size = 1) +
  geom_point(data = si_partition.df2, aes(y = fct_reorder(category, mean), x = value), size = 2, alpha = 0.7) +
  geom_line(data = si_partition.df2, aes(y = fct_reorder(category, mean), x = value, group = id), alpha = 0.2) +
  labs(x = "Fitness", y = "Conditions") +
  theme_bw()

Error in geom_vpline(data = si_partition.df2, aes(y = fct_reorder(category,  : 
  could not find function "geom_vpline"

#——- consequences of no targeted immunity ————# # get dynamics of no targeted immunity

get_dyn <- function(df){
  
  source(here("functions/chabaudi_si_clean_W.R"))
  id <- df$id
  cue <- df$cue
  log <- df$log
  par <- c(df$var_W1, df$var_W2, df$var_W3, df$var_W4)
  cue_range <- seq(df$low, df$high, by = df$by)
  
  # get dynamics
  dyn <- chabaudi_si_clean_W(
    parameters_cr = par,
    immunity = "tsukushi",
    parameters = parameters_tsukushi,
    time_range = seq(0, 20, by = 1e-3),
    cue_range =  cue_range,
    cue = cue,
    log_cue = log,
    solver = "vode",
    dyn = T
  )
  
  # combine
  dyn2 <- cbind(dyn, id = id, cue = cue, log = log)
  
  write_parquet(dyn2, paste0(here("data/partition/si_dyn/"), id, "_noW_dyn.parquet"))
  
}

get df to run

mclapply(si_partition.ls, get_dyn)

$`1`

$`2`

$`3`

$`4`

$`5`

$`6`

$`7`

$`8`

$`9`

$`10`

NA

process dataframe

# import in dataframe
no_W.ls <- list.files(here("data/partition/si_dyn/"), pattern = "*noW_dyn.parquet", full.names = T)
no_W.df <- lapply(no_W.ls, read_parquet)
no_W.df <- do.call(rbind, no_W.df)

# combine with ez label
ez_label <- read.csv(here("data/ez_label.csv"))
no_W.df <- left_join(no_W.df, ez_label, by = c("id" = "id_si"))

# get conversion rate 
no_W.cr <- no_W.df %>% filter(variable == "cr")
no_W.I <- no_W.df %>% filter(variable == "I")

# get default conversion rate dynamics
si_dyn.df <- left_join(si_dyn.df, ez_label, by = c("id" = "id_si"))
si_dyn.cr <- si_dyn.df %>% filter(variable == "cr")
si_dyn.I <- si_dyn.df %>% filter(variable == "I")

plot conversion rate

mechanism: targeted immunity led to lower parasite density in the initial stages, which prevents parasites from making the switch from no conversion rate to high conversion rate. When parsite density undergoes drastic increase at the beginning due to lower immunity, this presents a higher degree of signal that allows parasite to make the switch appropriately

partition_cr.pl <- ggplot() +
  geom_line(data = no_W.cr, aes(x = time, y= value, color = "No targeted immunity")) +
  geom_line(data = si_dyn.cr, aes(x = time, y= value, color = "Default")) +
  facet_wrap(~ez_label_si, ncol = 5) +
  xlim(0, 20) +
  geom_vline(xintercept = 7) +
  labs(x = "Time (days)", y = "Conversion rate", color = "Condition") +
  theme_bw()

no_W.cr

#—– cue state ————–#

function to get cue states

run function

plot

absence of targeted immunity led to drastic increase in parasite density in early phases of infection. This produces high signal intensity for parasite and host-based cues, especially non-logged ones, which allows for state differentation. While this can be viewed as a modelling artifiact, it should be noted that the logged cues seldom changed as these changes in early infection did little to alter the actual early signal intensity sensed by the parasite. In an environment where there is heterogeneity in host response, and thus, signal, logging allows for parasites to adapt optimal strategy whereas non-logged cues must contend with sensitivity to immunity.

  # arrange
  pl <- ggarrange(state_pl, cr_pl, ncol = 1, nrow = 2, align = "v", heights = c(1, 0.3))

Warning: Removed 4 rows containing missing values (geom_raster).

split

# combine state
noW_default.state <- left_join(
  select(no_W.state2, time, `No targeted\nimmunity` = value, id, ez_label_si), 
  select(default.state %>% filter(time <= 20), time, `Default` = value, id, ez_label_si), by = c("time", "id", "ez_label_si"))

noW_default.state2 <- tidyr::pivot_longer(noW_default.state, c(`No targeted\nimmunity`, `Default`))
# combine conversion raster
noW_default.cr <- left_join(
  select(no_W.cr, time, `No targeted\nimmunity` = value, id, ez_label_si), 
  select(si_dyn.cr %>% filter(time <= 20), time, `Default` = value, id, ez_label_si), by = c("time", "id", "ez_label_si"))
noW_default.cr2 <- tidyr::pivot_longer(noW_default.cr, c(`No targeted\nimmunity`, `Default`))

# split
noW_default_state.ls <- split(noW_default.state2, noW_default.state2$id)
noW_default_cr.ls <- split(noW_default.cr2, noW_default.cr2$id)

# run function
mapply(plot_state, noW_default_state.ls, noW_default_cr.ls)

#——– reaction norms of default vs optimized ————# # get reaction norm and rug data

source(here("functions/par_to_df.R"))

# Gametocyte
g_log.rn <- par_to_df(par = c(1.211521, -3.936778,  -1.312944,  -1.285713), cue_range = seq(0, log10(6*(10^4)), by = (log10(6*(10^4)))/5000))
g_log.rn2 <- par_to_df(par = c(1.393860539, -4.253007616,   -0.313947029,   -2.000857344), cue_range = seq(0, log10(6*(10^4)), by = (log10(6*(10^4)))/5000))

g.rn <- par_to_df(par = c(0.04061288,   -9.31445958,    74.13015506,    -431.5984364), cue_range = seq(0, 6*(10^4), by = (6*(10^4))/5000))
g.rn2 <- par_to_df(par = c(0.541729073, -3.904616443,   0.87487412, -0.694177021), cue_range = seq(0, 6*(10^4), by = (6*(10^4))/5000))

# I+Ig
I_Ig_log.rn <- par_to_df(par = c(3.594042,  4.157744,   -13.530672, 2.599905), cue_range = seq(0, log10(6*(10^6)), by = (log10(6*(10^6)))/5000))
I_Ig_log.rn2 <- par_to_df(par = c(63.71893822,  -87.77671601,   -56.55475514,   -66.02209549), cue_range = seq(0, log10(6*(10^6)), by = (log10(6*(10^6)))/5000))

I_Ig.rn <- par_to_df(par = c(0.3159297, -46.1104558,    1250.752908,    -6.1982093), cue_range = seq(0, 6*(10^6), by = (6*(10^6))/5000))
I_Ig.rn2 <- par_to_df(par = c(0.731982784,  -21.69799449,   149.7841876,    17.02551769), cue_range = seq(0, 6*(10^6), by = (6*(10^6))/5000))

# convert log to non-logged scale
g_log.rn$cue_range <- 10^(g_log.rn$cue_range)
g_log.rn2$cue_range <- 10^(g_log.rn2$cue_range)
I_Ig_log.rn$cue_range <- 10^(I_Ig_log.rn$cue_range)
I_Ig_log.rn2$cue_range <- 10^(I_Ig_log.rn2$cue_range)

# get rug
g_log.rug <- default.state %>% 
  filter(label_si == "G log") %>% 
  mutate(value = 10^value) %>% 
  select(label_si, value)

g_log.rug2 <- no_W.state %>% 
  filter(label_si == "G log") %>% 
  mutate(value = 10^value) %>% 
  filter(value <= 6*(10^4)) %>% 
  select(label_si, value)

I_Ig_log.rug <- default.state %>% 
  filter(label_si == "I+Ig log") %>% 
  select(label_si, value)

I_Ig_log.rug2 <- no_W.state %>% 
  filter(label_si == "I+Ig log") %>% 
  select(label_si, value)

g.rug <- default.state %>% 
  filter(label_si == "G") %>% 
  select(label_si, value)

g.rug2 <- no_W.state %>% 
  filter(label_si == "G" & value <= 6*(10^4)) %>% 
  select(label_si, value)

I_Ig.rug <- default.state %>% 
  filter(label_si == "I+Ig") %>% 
  select(label_si, value)

I_Ig.rug2 <- no_W.state %>% 
  filter(label_si == "I+Ig") %>% 
  select(label_si, value)

# get rug limits
rug_lim <- rbind(g_log.rug,
                 g_log.rug2,
                 I_Ig_log.rug,
                 I_Ig_log.rug2,
                 g.rug,
                 g.rug2,
                 I_Ig.rug,
                 I_Ig.rug2) %>% 
  group_by(label_si) %>% 
  summarize(max = max(hablar::s(value), na.rm = T),
            min = min(hablar::s(value), na.rm = T))

# combine and filter
rn <- rbind(
  cbind(g_log.rn, label_si = "G log", condition = "Default"),
  cbind(g_log.rn2, label_si = "G log", condition = "No targeted\nimmunity"),
  cbind(g.rn, label_si = "G", condition = "Default"),
  cbind(g.rn2, label_si = "G", condition = "No targeted\nimmunity"),
  cbind(I_Ig_log.rn, label_si = "I+Ig log", condition = "Default"),
  cbind(I_Ig_log.rn2, label_si = "I+Ig log", condition = "No targeted\nimmunity"),
  cbind(I_Ig.rn, label_si = "I+Ig", condition = "Default"),
  cbind(I_Ig.rn2, label_si = "I+Ig", condition = "No targeted\nimmunity")
) %>% 
  left_join(rug_lim, by = "label_si") %>% 
  group_by(label_si) %>% 
  filter(cue_range <= max & cue_range >= min)

# combine rug
rug <- rbind(cbind(g_log.rug, condition = "Default"),
             cbind(g_log.rug2, condition = "No targeted\nimmunity"),
             cbind(g.rug, condition = "Default"),
             cbind(g.rug2, condition = "No targeted\nimmunity"),
             cbind(I_Ig_log.rug, condition = "Default"),
             cbind(I_Ig_log.rug2, condition = "No targeted\nimmunity"),
             cbind(I_Ig.rug, condition = "Default"),
             cbind(I_Ig.rug2, condition = "No targeted\nimmunity"))

# cobine with ezlabel
rn2 <- rn %>% left_join(ez_label, by = "label_si")
rug2 <- rug %>% left_join(ez_label, by = "label_si")

# filter rug
default.rug <- rug2 %>% filter(condition == "Default")
no.rug <- rug2 %>% filter(condition == "No targeted\nimmunity")

plot

ggplot() +
  geom_line(data = rn2, aes(x = cue_range, y = cr, color = condition)) +
  geom_rug(data = default.rug, aes(x = value, color = condition), sides = "b") +
  geom_rug(data = no.rug, aes(x = value, color = condition), sides = "t") +
  facet_wrap(~fct_relevel(ez_label_si, c("Gametocyte log10", "Gametocyte", "Asexual&sexual\niRBC log10", "Asexual&sexual iRBC")), scales = "free_x") +
  scale_x_continuous(labels = function(x) format(x, scientific = TRUE)) +
  theme_bw() +
  scale_color_manual(values =c("Default" = "#4575b4", "No targeted\nimmunity" = "#fc8d59")) +
  ylim(0, 1.1) +
  labs(x = "Cue range", y = "Conversion rate", color = "Condition")

ggsave(here("figures/plos-bio/partition_rn.tiff"), width = 7.5, height = 6)

get conversion rate legend

noW_default.cr %>% filter(id == "G_log") %>% 
ggplot() +
  geom_raster(aes(x = time, y = id, fill = Default)) +
  xlim(1,20) +
  theme_bw() +
    labs(x = "Time (days)",
         fill = "Conversion rate") +
  theme(axis.title.y=element_blank(),
        axis.ticks.y=element_blank(),) +
    scale_fill_viridis_c(lim = c(0, 1))

ggsave(here("figures/plos-bio/cr_legend.tiff"))

#——————————————-# Best cue adoption consequences #——————————————-# # get data for disease curves

# single infection dynamics
si_dyn.df <- read_parquet(here("data/si_dyn/si_dyn_30.parquet")) 

# co-infection dynamics (mon-cue)
ci_dyn.df <- read_parquet(here("data/ci_dyn/ci_dyn.parquet"))

# dual cue dynamics
dual_dyn.df <- read_parquet(here("data/dual_cue_dyn/dual_cue_dyn.parquet"))

#——- single cue comparison —————# # process data

write_parquet(si_dc.low, here("data/disease_curve/si_dc_low.parquet"))

Error in write_parquet(si_dc.low, here("data/disease_curve/si_dc_low.parquet")) : 
  object 'si_dc.low' not found

plot

# plot
si_dc.pre <- ggplot() +
  geom_path(data = si_dc.poor, aes(x= total, y = R, group = id), color = "dark grey", arrow = arrow(type = "closed", angle = 10, length = unit(0.2, "inches"))) +
  labs(color = "Single infection\nhigh-performing cues", x = "Asexual & sexual iRBC", y = "RBC") +
  theme_bw()+ scale_y_continuous(labels = function(x) format(x, scientific = TRUE, accuracy = 0.1)) +
  guides(shape = FALSE)

Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

#———- co-infection monocue ————-#

# get relevent variables
ci_dc.df <- ci_dyn.df %>% 
  filter(variable == "I1" | variable == "Ig1" | variable == "R")

# morph into skinny format
ci_dc.df <- tidyr::pivot_wider(ci_dc.df, names_from = variable, values_from = value, id_cols = c(time, label)) %>% 
  mutate(total = I1+Ig1)

# good cue bad cue
ci_cue.dv <- ci_fitness.df %>% 
  mutate(classification = case_when(
    value > 2.25 ~ "High-performing",
    value <= 2.25 ~ "Poor-performing"
  ))

# join with classificaiton
ci_dc.df2 <- ci_dc.df %>% left_join(ci_cue.dv, by = "label")

# split into top erforming and poor-performing cues
ci_dc.high <- ci_dc.df2 %>% filter(classification == "High-performing")
ci_dc.poor <- ci_dc.df2 %>% filter(classification == "Poor-performing")

# join high performing with label
ci_dc.high2 <- ci_dc.high %>% left_join(ez_label, by = c("label" = "label_ci"))

Error in is.data.frame(y) : object 'ez_label' not found

plot

# plot
ci_dc.pre <- ggplot() +
  geom_path(data = ci_dc.poor, aes(x= total, y = R, group = label), color = "dark grey", arrow = arrow(type = "closed", angle = 10, length = unit(0.2, "inches"))) +
  labs(color = "Co-infection\nhigh-performing cues", x = "Asexual & sexual iRBC", y = "RBC") +
  theme_bw()+ scale_y_continuous(labels = function(x) format(x, scientific = TRUE, accuracy = 0.1)) +
  guides(shape = FALSE)

Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

#——— dual cue ————————–# # process data

# turn skinny
dual_dc.df <- dual_dyn.df %>% 
  mutate(label_alt = paste(label, "+" , label_b)) %>% 
  select(label_alt, time, variable, value) %>% 
  filter(variable == "I" | variable == "Ig" | variable == "R") %>% 
  distinct(label_alt, time, variable, .keep_all = T)

dual_dc.df2 <- dual_dc.df %>% 
  tidyr::pivot_wider(names_from = variable, values_from = value, id_cols = c(time, label_alt)) %>%
  mutate(total = I+Ig)

write_parquet(dual_dc.df2, here("data/disease_curve/dual_dc.parquet"))

# good dual cue -> list of good performing dual cues that encompass wide variety of cues
selected_dual_cue <- c("R log + I log", "R + Ig log", "G log + I log", "G log + Ig log", "Ig + I log")
bad_dual_cue <- c("G + I", "R + Ig", "R log + Ig", "G + R", "G + R log", "G + Ig", "Ig + I", "R + I", "R log + I")

# get classification -> R log10 + I log10 as the only good one
dual_dc.high <- dual_dc.df2 %>% filter(label_alt %in% selected_dual_cue) %>% 
  mutate(label_alt = gsub("log", "log10", label_alt))
dual_dc.poor <- dual_dc.df2 %>% filter(label_alt %in% bad_dual_cue) %>% 
  mutate(label_alt = gsub("log", "log10", label_alt))
#write_parquet(dual_dc.high, here("data/disease_curve/dual_dc_high.parquet"))
#write_parquet(dual_dc.poor, here("data/disease_curve/dual_dc_poor.parquet"))

plot

dual_dc.pre <- ggplot() +
  geom_path(data = dual_dc.poor, aes(x= total, y = R, group = label_alt), color = "dark grey", arrow = arrow(type = "closed", angle = 10, length = unit(0.2, "inches"))) +
  labs(color = "High-performing\ndual cues", x = "Asexual & sexual iRBC", y = "RBC") +
  theme_bw()+ scale_y_continuous(labels = function(x) format(x, scientific = TRUE, accuracy = 0.1)) +
  guides(shape = FALSE)

Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

#——— co-infection static —————-#

# import in dynamics data
static_dyn.ls <- list.files(path = here("data/ci_static/"), pattern = "*.parquet", full.names = T)
static_dyn.ls <- lapply(static_dyn.ls, read_parquet)

# filter variable and transform
static_dyn.ls2 <- mclapply(static_dyn.ls, function(x){
  x %>% 
    filter(variable == "I1" | variable == "Ig1" | variable == "I2" | variable == "Ig2" | variable == "R") %>% 
    mutate(id_alt = paste(id_1, id_2)) %>% 
    tidyr::pivot_wider(names_from = variable, values_from = value, id_cols = c(time, id_alt)) %>%
  mutate(total1 = I1+Ig1, total2 = I2+Ig2)
})

static_dc.df <- do.call(rbind, static_dyn.ls2)
static_dc.df <- static_dc.df %>% 
  mutate(id_1 = gsub(" .*", "", id_alt),
         id_2 = gsub(".* ", "", id_alt)) %>% 
  filter(id_1 != id_2)
#write_parquet(static_dc.df, here("data/disease_curve/static_dc.parquet"))

further processing

static_dc.df <- read_parquet(here("data/disease_curve/static_dc.parquet"))
# get winners and losers
## import in fitness
static_fitness.df <- read.csv(here("data/ci_static.csv"))
## get winner situation
static_fitness.df2 <- static_fitness.df %>% 
  filter(id_1 != id_2) %>% 
  mutate(winning_id = case_when(
    fitness_difference > 0 ~ id_1,
    fitness_difference< 0 ~ id_2
  ),
  losing_id = case_when(
    fitness_difference < 0 ~ id_1,
    fitness_difference> 0 ~ id_2
  ))

# left join
static_dc.df2 <- static_dc.df %>% 
  left_join(select(static_fitness.df2, id_1, id_2, winning_id, losing_id, fitness_difference), by = c("id_1", "id_2"))

# get winner-loser difference in terms of I+Ig also filter out to onyl very strong fitness difference
static_dc.df3 <- static_dc.df2 %>% 
  mutate(total_diff = case_when(
    fitness_difference > 0 ~ total1-total2,
    fitness_difference< 0 ~ total2-total2
  ),
  total_winner = case_when(
    fitness_difference > 0 ~ total1,
    fitness_difference< 0 ~ total2
  ),
  total_loser = case_when(
    fitness_difference > 0 ~ total2,
    fitness_difference< 0 ~ total1
  )) %>% 
  filter(abs(fitness_difference) > 0.5)

plot

static_dc.pl <- ggplot() +
  geom_path(data = static_dc.df3, aes(x= total_winner, y = R, group = id_alt, color = "Winner"), alpha = 0.5, arrow = arrow(type = "closed", angle = 10, length = unit(0.2, "inches"))) +
  geom_path(data = static_dc.df3, aes(x= total_loser, y = R, group = id_alt, color = "Loser"),
          alpha = 0.5,arrow = arrow(type = "closed", angle = 10, length = unit(0.2, "inches"))) +
  labs(color = "Status", x = "Asexual & sexual iRBC", y = "RBC") +
  scale_color_manual(values=c("Winner" = "#4575b4","Loser"= "#fc8d59"))  +
  theme_bw() + 
  scale_y_continuous(labels = function(x) format(x, scientific = TRUE, accuracy = 0.1))

#———co-infection invasion —————# # get invasion dynamic

# get invasion df
invasion_fitness.df <- read.csv(here("data/ci_invasion.csv"))

# get cue range
ci_cue_range <- read.csv(here("data/cue_range_ci.csv"))
invasion_fitness.df2 <- invasion_fitness.df %>% 
  left_join(select(ci_cue_range, id, mut_cue = cue, mut_low = low, mut_high = high, mut_by = by), by = c("mut_id"= "id")) %>% 
  left_join(select(ci_cue_range, id, res_cue = cue, res_low = low, res_high = high, res_by = by), by = c("res_id"= "id"))

function to get dynamic

get_invasion_dyn <- function(df){
  # get cues
  mut_cue <- df$mut_cue
  res_cue <- df$res_cue
  
  # get info of cues (for co infection)
  if(stringr::str_detect(mut_cue, "-i")){mut_cue = gsub("*-i", "1", mut_cue)}
  if(stringr::str_detect(mut_cue, "-i", negate = T)){mut_cue = mut_cue}
  if(stringr::str_detect(res_cue, "-i")){res_cue = gsub("*-i", "2", res_cue)}
  if(stringr::str_detect(res_cue, "-i", negate = T)){res_cue = res_cue}
  
  # get log
  mut_log <- ifelse(stringr::str_detect(df$mut_id, "log"), "log10", "none")
  res_log <- ifelse(stringr::str_detect(df$res_id, "log"), "log10", "none")
  
  # get parameters
  mut_par <- c(df$mut_var1_opt, df$mut_var2_opt, df$mut_var3_opt, df$mut_var4_opt)
  res_par <- c(df$res_var1, df$res_var2, df$res_var3, df$res_var4)
  
  # get cue range
  mut_cue_range <- seq(df$mut_low, df$mut_high, by = df$mut_by)
  res_cue_range <- seq(df$res_low, df$res_high, by = df$res_by)
  
  # get dynamics of co infection
  ci_dyn <- chabaudi_ci_clean(
    parameters_cr_1 = mut_par,
    parameters_cr_2 = res_par,
                  immunity = "tsukushi",
                  parameters = parameters_tsukushi,
                  cue_1 = mut_cue,
                  cue_2 = res_cue,
                  cue_range_1 = mut_cue_range,
                 cue_range_2 = res_cue_range,
                log_cue_1 = mut_log,
                log_cue_2 = res_log,
                solver = "vode",
                time_range = seq(0, 30, 0.001),
    dyn = T)
  
  # append label to all df
 ci_dyn2 <- cbind(ci_dyn, mut_id = df$mut_id, res_id = df$res_id)
  
  # write
 write_parquet(ci_dyn2, paste0(here("data/ci_invasion_dyn/"), df$mut_id, "-", df$res_id, ".parquet"))
}

run dynamic funciton

# get function and parameters
source(here("functions/chabaudi_ci_clean.R"))
parameters_tsukushi <- c(R1 = 8.89*(10^6), # slightly higher
                lambda = 3.7*(10^5),
                mu = 0.025, 
                p = 8*(10^-6), # doubled form original
                alpha = 1, 
                alphag = 2, 
                beta = 5.721, 
                mum = 48, 
                mug = 4, 
                I0 = 43.85965, 
                Ig0 = 0, 
                a = 150, 
                b = 100, 
                sp = 1,
                psin = 16.69234,
                psiw = 0.8431785,
                phin = 0.03520591, 
                phiw = 550.842,
                iota = 2.18*(10^6),
                rho = 0.2627156)
# split
invasion.ls <- split(invasion_fitness.df2, seq(nrow(invasion_fitness.df2)))

# run function
mclapply(invasion.ls, get_invasion_dyn, mc.cores = 4)

process data

# import in invasion dynamics
invasion_dyn.ls <- list.files(path = here("data/ci_invasion_dyn"), pattern = "*.parquet", full.names = T)
invasion_dyn.ls <- lapply(invasion_dyn.ls, read_parquet)

# filter and so on
invasion_dyn.ls2 <- mclapply(invasion_dyn.ls[167:177], mc.cores = 4, function(x){
  x2 <- x %>% 
    filter(variable == "I1" | variable == "Ig1" | variable == "I2" | variable == "Ig2" | variable == "R") %>% 
    mutate(id_alt = paste(mut_id, res_id)) %>% 
    select(id_alt, time, variable, value) %>% 
    tidyr::pivot_wider(names_from = variable, values_from = value, id_cols = c(time, id_alt)) %>%
  mutate(total1 = I1+Ig1, total2 = I2+Ig2)
  
  write_parquet(x2, paste0(here("data/disease_curve/ci_invasion/"), unique(x2$id_alt), "_dc.parquet"))
})

# fetch data
invasion_dyn.ls2 <- list.files(path = here("data/disease_curve/ci_invasion"), pattern = "*.parquet", full.names = T)
invasion_dyn.ls2 <- lapply(invasion_dyn.ls2, read_parquet)
invasion_dc.df <- do.call(rbind, invasion_dyn.ls2)
invasion_dc.df <- invasion_dc.df %>% 
  mutate(mut_id = gsub(" .*", "", id_alt),
         res_id = gsub(".* ", "", id_alt)) %>% 
  filter(mut_id != res_id)
#write_parquet(invasion_dc.df, here("data/disease_curve/invasion_dc.parquet"))

further processing

invasion_dc.df <- read_parquet(here("data/disease_curve/invasion_dc.parquet"))
# get winners and losers
invasion_fitness.df <- read.csv(here("data/ci_invasion.csv"))
invasion_dc.df2 <- invasion_dc.df %>% 
  left_join(invasion_fitness.df, by = c("mut_id", "res_id")) %>% 
  mutate(
  total_winner = case_when(
    fitness> 0 ~ total1,
    fitness< 0 ~ total2
  ),
  total_loser = case_when(
    fitness > 0 ~ total2,
    fitness < 0 ~ total1
  )) %>% 
  filter(abs(fitness) > 0.5)

plot

invasion_dc.pl <- ggplot() +
  geom_path(data = invasion_dc.df2, aes(x= total_winner, y = R, group = id_alt, color = "Winner"), alpha = 0.5, arrow = arrow(type = "closed", angle = 10, length = unit(0.2, "inches"))) +
  geom_path(data = invasion_dc.df2, aes(x= total_loser, y = R, group = id_alt, color = "Loser"),
          alpha = 0.5,arrow = arrow(type = "closed", angle = 10, length = unit(0.2, "inches"))) +
  labs(color = "Status", x = "Asexual & sexual iRBC", y = "RBC") +
  scale_color_manual(values=c("Winner" = "#4575b4","Loser"= "#fc8d59"))  +
  theme_bw() + 
  scale_y_continuous(labels = function(x) format(x, scientific = TRUE, accuracy = 0.1)) %>% 
  theme(legend.position = "none")

Error in `merge_element()`:
! No method for merging ScaleContinuousPosition into element_line
Backtrace:
 1. ggplot2:::`+.gg`(...)
 2. ggplot2:::add_ggplot(e1, e2, e2name)
 4. ggplot2:::ggplot_add.theme(object, p, objectname)
 5. ggplot2:::add_theme(plot$theme, object)
 7. ggplot2:::merge_element.default(t2[[item]], t1[[item]])

#——— plot disease curves together ————#

dc.pl1 <- ggarrange(si_dc.pl, ci_dc.pl, dual_dc.pl, ncol = 3, align = "h", labels = c("A", "B", "C"), widths = c(1.1, 0.9, 1.1))

Warning: stack imbalance in 'lapply', 97 then 98
Warning: stack imbalance in 'lapply', 84 then 85
Warning: stack imbalance in 'lapply', 71 then 72

#——— quantifying disease curve area ————# # function to calculate area between sets of points -> from https://stackoverflow.com/questions/3672260/area-covered-by-a-point-cloud-with-r

library(splancs)
cha<-function(df){
  x <- df$total
  y <- df$R
chull(x,y)->i
return(areapl(cbind(x[i],y[i])))
}

loop to get area: single infection

coinfection

# edit and join
ci_dc_high.area2 %>% ci_dc_high.area %>% 
  select(area, value) %>% 
  mutate(condition = "Co-infection")

Error in ci_dc_high.area(.) : could not find function "ci_dc_high.area"

dual cue

#—— get fitted scatter plot for all single infection, co infection, and dual cue ——–#

# plot
ggplot() +
  geom_point(data = dc.area, aes(x = area, y = value)) +
  facet_wrap(~condition, scales = "free") +
  geom_smooth(data = dc.area, method = "lm", alpha = .15, aes(x = area, y = value), color = "#4575b4") +
  labs(x = "Area", y = "Fitness", color = "Status") +
  theme_bw()

`geom_smooth()` using formula 'y ~ x'

#——- plot together with disease curve ——–#

ggarrange(si_dc.pl, ci_dc.pl, dual_dc.pl, dc_area.pl1, labels = c("A", "B", "C", "D"), align = "h", ncol = 4, widths = c(1, 0.9, 1, 1.5))

`geom_smooth()` using formula 'y ~ x'
Warning: Graphs cannot be horizontally aligned unless the axis parameter is set. Placing graphs unaligned.

ggarrange(si_dc.pl, ci_dc.pl, dual_dc.pl, dc_area.pl1, labels = c("A", "B", "C", "D"), align = "h", ncol = 4, widths = c(1, 0.9, 1, 1.5))

`geom_smooth()` using formula 'y ~ x'
Warning: Graphs cannot be horizontally aligned unless the axis parameter is set. Placing graphs unaligned.

ggsave(here("figures/plos-bio/disease_curve1.tiff"), units = "px", width = 2250, height = 600, scale = 2.5, dpi=300, compression = "lzw")

#——— static area comparison ————-# # compute area

# get winner and loser
static_dc.df4 <- static_dc.df %>% 
  left_join(select(static_fitness.df, id_1, id_2, fitness_difference), by = c("id_1", "id_2")) %>%
  filter(id_1 != id_2) %>% 
  mutate(
  total_winner = case_when(
    fitness_difference > 0 ~ total1,
    fitness_difference< 0 ~ total2
  ),
  total_loser = case_when(
    fitness_difference > 0 ~ total2,
    fitness_difference< 0 ~ total1
  ))%>% 
  filter(fitness > 0.5)

Error in `h()`:
! Problem with `filter()` input `..1`.
ℹ Input `..1` is `fitness > 0.5`.
x object 'fitness' not found
Backtrace:
 1. ... %>% filter(fitness > 0.5)
 7. base::.handleSimpleError(...)
 8. dplyr h(simpleError(msg, call))

plot static

#——— invasion area comparison —————–# # get area

# pair
invasion.area <- cbind(select(invasion_win.area, Winner = area),
                     select(invasion_loser.area, Loser = area)) %>% 
  mutate(classification = ifelse(Winner>Loser, "Winner larger area", "Loser larger area"))

Error in select(invasion_loser.area, Loser = area) : 
  object 'invasion_loser.area' not found

plot

#—— plot together ————-#